* Coal inventory regressions for "Green Expectations: Current Effects of Anticipated Carbon Pricing".

clear
capture log close
estimates clear

local dir="C:/yourdirectory"  /* main directory */
local dir_input = "`dir'"  /* directory with covariate files */

local dir_output="`dir'/output/coalstocks"  /* main directory */

log using "`dir'/logs/coalstocks.log", replace

set more off
set matsize 1000


*_______________________________________________________________________________


/* Get electricity sector and total coal stocks (thousand short tons)
   All values are for the end of the indicated period  				*/
use "`dir_input'/eia_coalstocks.dta", clear

tsset Counter

gen delta_electr = Electricity-L.Electricity
gen delta_total = Total-L.Total
gen pct_electr = ln(Electricity/L.Electricity)
gen pct_total = ln(Total/L.Total)

gen eventmonth = 0
replace eventmonth = 1 if (Year==2010 & Month==4)

* Electricity sector's stock
xi: ivreg2 pct_electr eventmonth, kernel(bartlett) bw(auto) small	
estimates store nofe
xi: ivreg2 pct_electr eventmonth i.Month, kernel(bartlett) bw(auto) small	
estimates store monthfe
xi: ivreg2 pct_electr eventmonth i.Month i.Year, kernel(bartlett) bw(auto) small	
estimates store yearfe


* Total U.S. stock - Not doing this because data other than electric power sector are imputed from quarterly data by using 1/3 of the quarterly change
/*
xi: ivreg2 pct_total eventmonth, kernel(bartlett) bw(auto) small	
xi: ivreg2 pct_total eventmonth i.Month, kernel(bartlett) bw(auto) small	
xi: ivreg2 pct_total eventmonth i.Month i.Year, kernel(bartlett) bw(auto) small	
*/

local specif_list "nofe monthfe yearfe"

foreach specif of local specif_list{
	/* obtain residuals so can do SQ test and normality tests */
	estimates restore `specif'
	predict resid_`specif', residual
	*egen std_resid_`specif'  = std(resid_`specif')

	* save residuals for taking into Excel
	save "`dir_output'/temp_results.dta", replace
	drop if eventmonth==1
	keep resid_*
	*xpose, clear varname
	save "`dir_output'/resids.dta", replace
	use "`dir_output'/temp_results.dta"  /* get all variables back */	
	
	drop if eventmonth==1
	* Shapiro-Wilk test for normality
	swilk resid_`specif'
	* Jarque-Bera test for normality
	jb6 resid_`specif'
	use "`dir_output'/temp_results.dta", clear  /* get all variables back */	
}

drop if eventmonth==1


/* Plot density of residuals and original variable*/
kdensity resid_yearfe, generate(x_resid_base y_resid_base) nograph			
keep x_* y_*
save "`dir_output'/kdensity_base.dta", replace
use "`dir_output'/temp_results.dta"  /* get all variables back */	
drop if eventmonth==1
kdensity pct_electr, generate(x_var y_var) nograph			
keep x_* y_*
save "`dir_output'/kdensity_originalvar.dta", replace
use "`dir_output'/temp_results.dta"  /* get all variables back */	

log close
